Gully erosion is an indicator of land degradation, which occurs due to the removal of soil along drainage channels through surface water runoff. Given the significant impact of gullies, it is essential to develop reliable and targeted analysis methods to better understand their spatial occurrence and enable improved representation of these erosion processes in catchment sediment budget models.
The aims of the study were:
Mangatu catchment for training and testing; Wairoa for model validation (Fig. 2.1).
Figure 2.1: Study areas for analysis. The Mangatu catchment was used for model training and testing and subsets of the Wairoa catchment for validation
The reference data consists of a combination of all gully features for 1939, 1957, 1960, 1970, 1988, 1997 from Marden et al. (2012) and Marden, Herzig, and Basher (2014) (Fig. 2.2). These data was manually delineated based on historical aerial photography.
Figure 2.2: Gully reference data for the Mangatu catchment. From Marden et al. (2012) and Marden, Herzig, and Basher (2014). Background: Aerial imagery collected between 2012 and 2013.
The original data delivered by MWLR had a problem with co-registration, where polygons were not matching their correct geographical locations. To correct this, we used an affine transformation, see more details here.
The data was split into training and testing with a 70:30 ratio, stratified by year. To use it as reference data, the polygons were merged and dissolved.
library(sf)
library(tidyverse)
library(here)
g = read_sf(here("data/reference/gullies.shp"))
set.seed(7319)
g = g %>% mutate(id = row_number())
train = g %>%
group_by(year) %>%
slice_sample(prop = 0.7) %>%
ungroup()
test = g %>%
as_tibble() %>%
anti_join(train, by = 'id') %>%
st_as_sf() %>%
select(-id)
train = train %>%
select(-id)
write_sf(
train,
dsn = here("data/reference/gullies_train.shp"),
delete_layer = TRUE
)
write_sf(
test,
dsn = here("data/reference/gullies_test.shp"),
delete_layer = TRUE
)
train %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf(class = 1) %>%
write_sf(
dsn = here("data/reference/gullies_train_union.shp"),
delete_layer = TRUE
)
test %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf(class = 1) %>%
write_sf(
dsn = here("data/reference/gullies_test_union.shp"),
delete_layer = TRUE
)
The LiDAR DEM data corresponds to the year 2019 with a spatial resolution of 1 m. We used SAGA to calculate 12 terrain derivatives based on the LiDAR DEM, which included:
We tested a region-based convolutional neural network (Mask-RCNN) deep learning approach for object detection to map gully features. Mask-RCNN is a model for image segmentation, developed on top of Faster R-CNN (He et al. 2017). See more information on how Mask-RCNN works here.
The deep learning was performed woth a RGB combination of the 12 terrain derivatives. Latest tests were acceptable with: ‘lsfct,’‘shade,’‘tridx.’ Labelled chips for training data were generated with the reference gully features.
The main workflow took place in ArcGIS Pro 2.9.2. The workflow in figure 2.3 is recommended. The operating system used was Microsoft Windows 10 Enterprise, with processor characteristics Intel(R) Xeon(R) CPU E5-1650 v4 @ 3.60GHz, 6 cores and 12 logical processors, 64GB of installed RAM and NVIDIA GeForce GTX 1070 GPU with 8GB of dedicated memory and 6.1 compute capability.
Figure 2.3: Deep Learning workflow for object detection in ArcGIS Pro
These steps can be applied using the GUI or with the ArcPy Python modules for ArcGIS.
Using the ArcPy and learn modules, we established the workflow below within a Jupyter notebook and a ArcGIS Pro Project.
See documentation of the Python API here.
See the whole Jupyter notebook here.
# Import system modules and check out ArcGIS Image Analyst extension license
import arcpy
arcpy.CheckOutExtension("ImageAnalyst")
from arcpy.ia import *
from arcpy import env
data_dir = r'data'
env.workspace = data_dir
# Composite bands
arcpy.CompositeBands_management(
"mangatu_cplan.tif;mangatu_cprof.tif;mangatu_ddgrd.tif;mangatu_lsfct.tif;mangatu_mbidx.tif;mangatu_shade.tif;mangatu_slope.tif;mangatu_spidx.tif;mangatu_textu.tif;mangatu_tridx.tif;mangatu_twidx.tif;mangatu_vdcnw.tif",
"terrain.tif"
)
# Rename bands
terrain = r'data\terrain.tif'
terrain_raster = Raster(terrain, True)
terrain_raster.renameBand(1, 'cplan')
terrain_raster.renameBand(2, 'cprof')
terrain_raster.renameBand(3, 'ddgrd')
terrain_raster.renameBand(4, 'lsfct')
terrain_raster.renameBand(5, 'mbidx')
terrain_raster.renameBand(6, 'shade')
terrain_raster.renameBand(7, 'slope')
terrain_raster.renameBand(8, 'spidx')
terrain_raster.renameBand(9, 'textu')
terrain_raster.renameBand(10, 'tridx')
terrain_raster.renameBand(11, 'twidx')
terrain_raster.renameBand(12, 'vdcnw')
terrain_raster.getRasterBands()
# Define variables
terrain = r'data\terrain.tif'
terrain_raster = Raster(terrain, True)
gullies = r'data\gullies_train_union.shp'
# Select bands
bands_sel = ['cplan','lsfct','slope','tridx','twidx']
terrain_sel = ExtractBand(terrain_raster, band_names=bands_sel)
terrain_sel.save('data/temp/terrain_'+''.join(bands_sel)+'.tif')
# Convert data to 8BIT
nd_val = -99999
terrain_8bit = "temp/terrain_" + ''.join(bands_sel) + "_8bit.tif"
arcpy.CopyRaster_management(
in_raster=terrain_sel,
out_rasterdataset=terrain_8bit,
nodata_value=nd_val,
pixel_type="8_BIT_UNSIGNED",
# For some reason, actually activating this option DOES NOT scale the pixels
scale_pixel_value="NONE",
format= "TIFF")
## Set dir for chip generation named after band selection
out_chips = 'work/terrain_'+ ''.join(bands_sel)
# Create output_dir if not existing
if not os.path.exists(out_chips):
os.mkdir(out_chips)
# Execute ArcPy tool for data export
ExportTrainingDataForDeepLearning(
terrain_sel, out_chips, gullies,
"TIFF", 1024, 1024, 512, 512,
"ONLY_TILES_WITH_FEATURES", "RCNN_Masks", 0, "",
0, None, 0, "MAP_SPACE", "PROCESS_AS_MOSAICKED_IMAGE",
"NO_BLACKEN", "FIXED_SIZE"
)
# Prepare data
from arcgis.learn import prepare_data
data = prepare_data(
out_chips,
dataset_type='RCNN_Masks',
batch_size=4, chip_size=512,
imagery_type='ms',
norm_pct=1
)
print(data)
data.show_batch(rgb_bands=[0,1,2])
# Define model architecture
from arcgis.learn import MaskRCNN
maskrcnn = MaskRCNN(
data, backbone='resnet50', pointrend=True
)
# Find efficient learning rate
maskrcnn.lr_find() # copy into next command
# Fit the model
maskrcnn.fit(epochs=20, ,early_stopping=True, lr=maskrcnn.lr_find())
# Unfreeze the model
maskrcnn.unfreeze()
# Plot losses
maskrcnn.plot_losses()
# Show results
maskrcnn.show_results(rows=3)
# Save model
model_name = '20ep'
maskrcnn.save(model_name)
Sadly, this section here does not work with code, so for now, please go to the Geoprocessing tool :)
arcpy.env.processorType = "GPU"
# Set local variables
terrain_sel = ExtractBand(terrain, band_names=['cplan','lsfct','slope','spidx','tridx','twidx'])
in_raster = terrain_sel
out_detected_objects = "work/terrain_3lay/results/gully_lsfctshadetridx_20ep.shp"
in_model_definition = "work/terrain_3lay/models/gully_model_lsfctshadetridx_20ep/gully_model_lsfctshadetridx_20ep.emd"
model_arguments = "padding 256;threshold 0.5;batch_size 4"
run_nms = "NO_NMS"
confidence_score_field = "Confidence"
class_value_field = "Class"
max_overlap_ratio = 0
processing_mode = "PROCESS_AS_MOSAICKED_IMAGE"
# Check out the ArcGIS Image Analyst extension license
arcpy.CheckOutExtension("ImageAnalyst")
# Execute
DetectObjectsUsingDeepLearning( in_raster, out_detected_objects,
in_model_definition, model_arguments, run_nms, confidence_score_field,
class_value_field, max_overlap_ratio, processing_mode)
For this step we tested several padding sizes (256, 128, 64, 32). Every run of the tool resulted in slightly different results.
To validate the deep learning model, gully features from the Wairoa catchment, northern Hawke’s Bay would be ideally used. This data is from the SedNetNZ project, where erosion features including gullies, earthflows and cliffs were mapped on a 5x5 km grid. The features were mapped based on aerial imagery. For the Northern Wairoa catchment, the Gisborne 2017-19 regional aerial imagery with 0.3 m spatial resolution and for the Southern Wairoa catchment, the HBRC 2019-20 regional aerial imagery with 0.3 m spatial resolution were used.
The validation would take place by selecting the mapped gully features from the dataset, and using them as a base for comparison with the results of the deep learning model applied for this area. Measures such as Intersection over Union (IoU) and accuracy can be then computed to evaluate the performance of the model using an unbiased and completely new dataset.
For now, work has been delayed due to the delivery of LiDAR for the Hawke’s Bay region (expected now in July 2022). However, some other LiDAR tiles were available covering two grid cells (62 and 119) from the erosion feature mapping. An initial exploration of these tiles are shown in figure 2.4.
Click to expand!